function res = iodine_reaction(t,V)

I = V(1);
S2O8 = V(2);
I2 = V(3);
SO4 = V(4);
S2O3 = V(5);
S4O6 = V(6);

k1 = 0.0039832/(0.2784^2);
k2 = 3000;

% 2I + S2O8 -> I2 + 2SO4
% I2 + 2S2O3 -> S4O6 + 2I

r1 = k1 * I^2 * S2O8;
r2 = k2 * S2O3 * I2;

dIdt = -r1 + r2;
dS2O8dt = -r1;
dI2dt = r1 - r2;
dSO4dt = r1;
dS2O3dt = -r2;
dS4O6dt = r2;

res = [dIdt ; dS2O8dt ; dI2dt ; dSO4dt ; dS2O3dt ; dS4O6dt];

end
